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Abstract 

We propose a method to construct numerical solutions of parabolic 
equations on the unit sphere. The time discretization uses Laplace 
transforms and quadrature. The spatial approximation of the solution 
employs radial basis functions restricted to the sphere. The method 
allows us to construct high accuracy numerical solutions in parallel. 
We establish L_2 error estimates for smooth and nonsmooth initial data, 
and describe some numerical experiments. 
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1 Introduction 

We consider the initial-value problem 

3 t u + Au = f(t), for t > 0, withu(0)=u , (1.1) 
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where 3 t = 3/3t and A is a linear, self-adjoint, positive-semidefinite, second- 
order elliptic partial differential operator on the unit sphere. In our standard 
example, —A is the Laplace-Beltrami operator. The source term f (t) may 
depend on the spatial variables but we suppress this dependence in our no- 
tation, viewing f (t) as an element of a function space on the sphere. 

Instead of using time stepping for the numerical solution, as was done 
previously |4j, our approach is to represent the solution of (II. ip as an inverse 
Laplace transform, which is then approximated by quadrature. Developed 
first for parabolic problems by Sheen, Sloan and Thomee such an ap- 
proach is also effective for some evolution equations with memory [5j . These 
and related papers have discussed thoroughly the time discretization, but 
for the space discretization have considered only piecewise linear finite el- 
ements on a bounded domain in K n . Here, we propose instead a space 
discretization using spherical radial basis functions (SRBFs), which are con- 
venient for parabolic problems on Riemannian surfaces such as the unit 
sphere § n = { x G M n+1 : |x| = 1 }. 

Denoting the Laplace transform of u with respect to t by 



u(z) =£{u{t)}:-- 



e~ zt u(t)dt, (1.2) 



o 



we find that the solution of fll.ip formally satisfies 

(zl + A)u(z) = g(z) := u + f (z), (1.3) 

where I denotes the identity operator. The spectrum of A is a subset of the 
half-line [0, oo), so if z ^ (— oo,0] and if the Laplace transform f(z) exists, 
then 

u(z) = (zI + A)- 1 g(z). (1.4) 

When f (z) is analytic and bounded for d\z > 0, the solution u(t) can be 
recovered via the Laplace inversion formula 



1 

ut = 

v ' 27ti 



e zt u(z)dz, fort>0, (1.5) 

r 



where To is the contour £Kz = w, for any w > 0, with Jz increasing. 

Section [2] summarizes some technical results and assumptions needed for 
our subsequent analysis. In Section [3] we describe the time discretization 
and quote a known error estimate (Theorem 13. ip . after which we introduce 
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the space discretization using SRBFs. The heart of the paper is Section HI 
where we prove two error bounds for the space discretization by adapting 
the analysis of Thomee |13| for a finite element approximation of the heat 
equation on a domain in IR n . The first bound (Theorem 14.51) requires some 
spatial regularity of Uo and f , and is proved by estimating a contour integral. 
The second bound is proved by an energy argument, and assumes f = but 
allows nonsmooth initial data u G L 2 (S n ). Both bounds include a factor 
that blows up as t 0. Finally, Section [5] describes the results of some 
numerical experiments. 



2 Preliminaries 

2.1 Resolvent estimates 

We now view A as an abstract, densely defined, self-adjoint and positive- 
semidefinite linear operator on a complex Hilbert space "K. Assume further 
that (1 + A) -1 : 0~C — > "K is compact, so A has a discrete spectrum, and order 
the eigenvalues ^ Ai ^ A 2 ^ • • • . Note that Aj — > oo as j — > oo if "K is 
infinite dimensional. 

For any cp > 0, the spectrum of A is a subset of a closed sector in the 
complex plane C, 

Z v :={z^ : |argz| ^ cp}U{0}, with < cp < n/2. 

In addition, there is a constant C > such that A satisfies the resolvent 
estimate || (zl — A) _1 || ^ C|z| _1 for z e C \ L v , or, equivalently, 

IKzI + A)- 1 !!^ C|z|-\ forzGl^^, (2.1) 
where || ■ || denotes the operator norm induced by the norm in "K. 



2.2 Sobolev spaces on the unit sphere 

Denote the inner product in !K = L.2(§ n ) by 



(v,w> := 



vw dS, 



where dS is the surface measure on the unit sphere, and denote the mea- 
sure of the whole sphere by cu n (so, for example, uo 2 = 47t). Recall jH] 
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that a spherical harmonic is the restriction to S n of a homogeneous polyno- 
mial Y(x) in M n+1 satisfying AY(x) = 0, where A is the Laplacian operator 
in M n+1 . The space of spherical harmonics of degree £, denoted by Jit, has 
dimension N(n, I) := dimJC^, given by 

N(n,0) = l and N(n. t) = (2< + n ~ t" = &r t > 1. 

£!(u — lj! 

In the usual way, we construct an orthonormal basis { Y Ck : 1 < k ^ N(n, £) } 
for JC £ , so that (Yf k , Yvw) = S««'6kk'- 

The Laplace-Beltrami operator A* on S n may be defined in terms of the 
Laplacian A on IR n+1 by 

A*v = Av|§n where v(x) = v(x/|x|). (2.2) 

The spherical harmonics are eigenfunctions of A*, satisfying 

-A*Y £k = A £ Y £k where A £ =£(£ + n-l), 

for 1 ^ k < N(n, I) and £ G {0, 1,2, . . .}. Every function v G L 2 (S n ) can be 
expanded in a generalized Fourier series 

oo N(u,«) 

v = ^ ^ V£ k Y« k where Vf k = (v,Y fk ), 

£=0 lc=l 

and for cr G K. we can characterize the Sobolev space on the unit sphere, 
H a = H a (S n ), in terms of the generalized Fourier coefficients: v G H a if and 
only if the norm defined by 

oo N(n,«) 

||v|| 2 m :=||(I-AT/ 2 v|| 2 = }jl + A £ r Y. l^l 2 ( 2 - 3 ) 

£=0 k=l 

is finite. We also define the subspace of functions with mean zero, 
H C =H T) :={vGH ff (§ n ): vdS = 0}; 

since Y 01 = l/^/cu n is constant, we see that v G H a belongs to Hq if and 
only if v i = 0. 
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2.3 Positive definite kernels on the unit sphere 

A continuous function O : § n x S n — y R is called a positive definite kernel |10[ 
ITT] on § n if it satisfies the following two conditions: 

(i) 0(x,y) = 0(y,x) for all x, y G § n ; 

(ii) for any set of distinct scattered points {yi,y2, • • • , yi<} C S n , the sym- 
metric K x K matrix [®(yt,yj)] is positive semi-definite. 

We call G> strictly positive definite if the matrix is strictly positive definite. 

We will work with a kernel G> defined in terms of a univariate function 4> : 
[-1, 1] R by 

0>(x,y) = cb(x--y) for all x, y G § n , (2.4) 

where x • y denotes the Euclidean inner product of x and y. Following 
Miiller [6], let P £ (t) denote the Legendre polynomial of degree I for R n+1 , 
and expand 4>(t) in a Fourier-Legendre series 

1 x ' 

*(t) = — y"N(n,£)a £ P £ (t). (2.5) 
Due to the addition formula for spherical harmonics [6, Page 10], 

Y_ Y £k (x)Y £k (y) = ^^ip,(x-y), 
the kernel O can be represented as 



N (Tt, 



<l>[x,y)=Y_ Y. a ^WY £k (y), (2.6) 



£=0 k=l 

and since P £ (l) = 1 we find that 



oo 

(D(x,-)||h- = ^(l + A £ ) T aj;N(n,£), for all x G § n . (2.7) 



Chen et al. [2J proved that the kernel O is strictly positive definite if and 
only if d£ ^ for all I ^ and a £ > for infinitely many even values of I 
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and infinitely many odd values of t; see also Schoenberg [ID] and Xu and 
Cheney [17J. Here, we assume there is a T > n/2 and positive constants 
c and C such that 

c(l + A£r T ^ a £ < C(1 + A £ )" T , for aU f > 0. (2.8) 

Hence, O is strictly positive definite and, since N(n, £) = 0(£ n_1 ) as £ — > oo, 
the sum (12.71) is finite so, for each fixed x £ S'" - , the function y i— >■ 0(x,xj) 
belongs to H T (S n ). Moreover, this function is continuous by the Sobolev 
imbedding theorem. 



3 The discrete problem 

Choose an angle (3 £ (n/2,n— cp) and let V be any curve in the interior 
of the sector Zp which is homotopic to the line T appearing in the Laplace 
inversion formula (ll.5p . Deforming the contour of integration in (jl.5j) . we 
may then write 

1 



ut = 



2m 



e zt u(z)dz, 



(3.1) 



assuming that f(z) is analytic on and to the right of T. 

By taking f = in (II. ip . so that g(z) = u in (II. 3p . we see that the solu- 
tion operator for the homogeneous problem has the integral representation 



2m 



e zt £(z)u dz, where £(z) = (zl + A) 



-i 



(3.2) 



For the inhomogeneous case, the inverse Laplace transform of £(z)f(z) is the 
convolution of £(t) and f(t), giving the Duhamel formula 



ut 



£(t)iio 



£(t-s)f(s)ds. 



(3.3) 



A standard energy argument shows that ||£(t)u || ^ ||uo|| f° r all t ^ 0, so 
the continuous problem (11.11) is stable in the sense that 



ft 



|u(t)|| ^ llttol 



||f(s)||ds, fort^O. 
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For our numerical methods we choose V to be the curve with parametric 
representation 

z{l) :=cu + A(l-sin(6-i£,)), for £, e R, (3.4) 
where the constants w, A and 5 satisfy 

w>0, A>0 and 0<5<(3 - n/2. (3.5) 
Writing z = x + ry, we find that F is the left branch of the hyperbola 

*-°>-^ f y ^ -1 (3.6) 



A sin 5 / \ A cos 6 

which cuts the real axis at the point z = a> + A(l — sin 5) and has asymptotes 
y = ±(x — cv — A) cot 5. Thus, the conditions (I3.5P ensure that V lies in the 
sector := tv + Zp C Ip, and crosses into the left half-plane. 

We use (13. 41) in ( 13. ip to represent u(t) as an integral with respect to 



u(t) = — 



w(z(£,))z'(£,)d£,. (3.7) 



Since |e z(t)t | = e^ 2 ^^ = e^V^ 1 - 81 " 5 ™ 811 ^, the integrand exhibits a double 
exponential decay as |£J — > oo, for any fixed t > 0. 

3.1 Time discretization 

We choose a quadrature step k, put 

and apply an equal weight rule to the integral (13. 7p to obtain an approximate 
solution 

k N 

U N (t):=— Y_ e Zit u(zj)z;. (3.8) 

j=-N 

In view of (13. ip . to compute Un (t) we must solve the 2N + 1 equations 

(z j I + A)u(z j ) = g(z j ), for|j|<N. (3.9) 

These equations are independent and hence may be solved in parallel. Notice 
that the u(zj) determine the approximate solution (13. 8p for all t > and that 
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the numerical solution (13. 8p depends on the choice of the curve P, even though 
the representation (13. ip does not. However, we will see that a given V and 
k yield an accurate approximation Un (t) ~ u(t) only for t at a particular 
time scale. 

The parametric representation (13 .4p of V extends to a conformal mapping 
z = W{Q = cu + A(l-sin(6-iO), (3.10) 

which, for r > 0, transforms the strip Y r := {C : ^ t} onto the set S T := 
{^(C) : C G Y r } D P. In fact, W maps the line 3C = t\ to the left branch of a 
hyperbola given by (13. 6p with 6 replaced by 5 +T|. Thus, S r is bounded by 
the left branches of the hyperbolas corresponding to 3£ = r and J£ = — r. 
To ensure that S r C Ljj and that fRz — > — oo if ]z| — > oo with z G S r , we 
require < 6 — r<6 + r< (3 — 7t/2, or equivalently that 

< r < min(6, (3 -tt/2-6). (3.11) 

We introduce the notation 

||g||x,z '■— sup ||g(z) ||x, for X C IK and Z C C, 

zez 

abbreviated by ||g||z if X = IK, and put lg(s) = max(l, log(l/s)) . 

Theorem 3.1. Let u be the solution of (II. ip . f bounded and analytic 
in and /ix a time sca/e T > 0. Let < 6 < 1 and define b > 
by coshb = 4/(9 sin 6), let r satisfy (13.1 ip so that Y C S r C Z^, and put 
A = 7tr9N/(bT). TTien the approximate solution UnU) defined by (13 .8p 
with k = b/N ^ 27trlog2 satisfies 

||U N (t) -u(t)|| < Ce^lgfp.Nje-^dKH + ||f|| £? ), /or T/2 ^ t < 2T, 

where \i = 27tr(l — 9)/b, p r = 7tr9 sin(5 — r)/(2b) and C = Cs^p. 

Proof. See McLean and Thomee Theorem 3.1]. □ 

3.2 Galerkin approximation by SRBFs 

Given a suitable set of points X = {xi,x%, . . . , xk} C S n and a strictly positive 
definite kernel ® (x, y ) , we define the spherical radial basis functions O p (x) := 
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0(Xp,x) for 1 ^ p ^ K. Recall that our assumption ( 12. 8 p ensures O p G H T 
with T > n/2 ^ 1; thus 

S h := span{O p : 1 <p ^ K}C H 1 . 

The uniformity of the set X is measured by its mesh norm fix and its sepa- 
ration radius qx, defined by 

H = fix := sup mincos -1 ^ • x) and q = q x := - min cos _1 ("y • x). 

In words, fix is the maximum geodesic distance from a point on § n to the 
nearest point of X. For our convergence analysis, we require that the family 
of point sets {X} has a bounded mesh ratio: 

H x < Cq x . (3.12) 

Associated with the second-order, partial differential differential opera- 
tor A is a bounded sesquilinear form a : H 1 x H 1 — > C defined by 

a(u,v) = (Au,v) for u, v G H 1 . 

For example, if A = —A* then a(u, v) = (gradu, gradv) where grad is the 
surface gradient. The mild solution u : [0, oo) — > L2(§ n ) of (II. ip satisfies 

(3 t u,v) + a(u, v) = (f(t),v) for t > and all v G H 1 , 

with u(0) = Uo, and we define a semidiscrete solution Uh : [0, oo) — > Sh 
of ALU) by 

(3 t u H ,x) + a(u H ,x) = <f(t),x) forallxeS H , (3.13) 

with Uh(0) = Uoh ~ no for a suitable Uoh G Sh- 

The Laplace transform of u at Zj is the weak solution u(zj ) G H 1 of ( 13. 9p . 
that is, 

zj(u(zj),v) + a(u(Zj),v) = (g(zj),v) for all v g H 1 , 
and the Laplace transform of the semidiscrete solution, Un(zj) G Sh, satisfies 
Zj(fth(zj),X> + a(u H (Zj),x) = (9h(zj),x) for all x e S H , (3.14) 
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where gh[z) = Uoh + Ph.f(z) G Sh and Ph denotes the orthogonal projector 
from L 2 (S n ) onto Sh- Thus, we can view Uh_(zj) as a Galerkin approximation 
to u(zj ). Concretely, to compute u h (z) = Y.p=i Up we form the KxK 

matrices B and S, with entries 

B pq = (<D p ,O q ) and S pq = a(<D p , O q ), (3.15) 

form the load vector G(z) G C K with components G p (z) = (gn(z), ® p ), and 
then solve the KxK complex linear system 

(z j B + S)U(z j ) = G(z j ), (3.16) 

to obtain the solution vector \l[z) G C K with components U p (z). In contrast 
to finite element mass and stiffness matrices, B and S are not sparse because 
the SRBFs have large supports. 

3.3 Fully-discrete solution 

Combining the time and space discretizations, we arrive at a fully-discrete 
solution 

k N 

Un.hW = -j; e^UHtzjJzj, (3.17) 

1 J=-N 

whose evaluation requires that we solve the linear system (I3.16P at each of 
the 21M + 1 quadrature points Zj . (In practice, we also use quadratures for the 
integrations over §> n that are needed to compute B pq , S pq and G p (z), but 
for our analysis we assume that these quantities are computed exactly.) The 
elliptic differential operator A induces a discrete operator : Sh — > Sh, 
defined by 

(Ah4>,x) = a(i|),x), for ip, x e S H , (3.18) 
and the Galerkin equations (I3.14p are equivalent to 

( Zj I + AH)uH(zj) = g h (zj). (3.19) 

If we choose Uoh = PhU-o then gh(zj) = Png(zj) and by taking CK = Sh 
equipped with the L 2 -norm, we can apply Theorem 13.11 to Ah and deduce 
that 

||U N ,h(t)-Uh(t)|| < Ce^lgfprNJe-^dluoll+HfH^), for T/2 ^ t < 21. 

(3.20) 
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Since the triangle inequality gives 

||U Njh (t) -u(t)|| < ||U N>H (t) -u H (t)|| + ||u H (t) -u(t)||, (3.21) 

to estimate the error in U^ k it now suffices to estimate the error in the 
semidiscrete approximation Un(t). 

4 Error analysis of the spatial discretization 

We assume now that A = —A*. Since A = but A^ ^ A x = n for all I ^ 1, 
we see that 1 + A^ ^ (1 + n _1 )A^ for all I ^ 1. Hence, the sesquilinear form a 
is coercive on Hq, that is, 

vdS = 0. (4.1) 

Our analysis follows Thomee [TBI Chapter 3], with A* in place of the Lapla- 
cian (with homogeneous Dirichlet boundary conditions). Some technical 
modifications are needed, however, because A* has a zero eigenvalue. 

4.1 Approximation by SRBFs 

We will use the following estimate for the best approximation by SRBFs. 

Theorem 4.1. Assume that the Fourier-Leg endre coefficients in the expan- 
sion (12. 5p satisfy (12. 8 j) with x > n/2, so that S h C H T (§ n ). For any real 
q and v satisfying q ^ v $C 2t and q ^ r, ifvE H v then there exists x £ Sh 
such that 

||X — v||hi < Ch£ _q ||v|| H v. 

Proof. See Tran et al. [T5| Theorem 3.2] or [141 Theorem 3.7 and Remark 5.1], 
and note our assumption (I3.12p . □ 

In the special case q = 0, the estimate must hold for x = Ph.v, giving the 
following result. 

Corollary 4.2. The ^-projection ofv onto Sh has the approximation prop- 
erty 

||v-P H v|| ^ CH£||v|| H v for0^y^2r. 



a(v, v) ^ — - if v G H 1 and v 10 = 

1 + 
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For our error analysis, we also use the Ritz projector R^ : H 1 (§ n ) — > 
determined by the sesquilinear form 

ai(u, v) = a(u,v) + (u,v) for u, v G H 1 . 

We see from (14. ip that ai is coercive on H 1 ; in fact, ai(v, v) = ||v||^i. Thus, 
RhV G Sh is well-defined by 

ai (R H v, x) = ai (v, x) for all x G S h , (4.2) 

and the following error estimates hold using standard arguments. 

Theorem 4.3. IfvE H v and 1 ^ v < 2t ; then 

||v — R K v|| H i = inf ||v — xIIhi < Ch.^ _1 ||v|| H v 
xes h 

and 

II v- R h v || ^ CHx||v||m- 
Proof. The definition (14. 2D immediately implies the orthogonality property 

a x (v - R h v, x) = for all x e S H , (4.3) 

so, because ai(v, v) = HvU^i, 

||v-R h v||hi = ai(v-R H v,v-R h v) = a x (v - R H v, v -x) 
< ||v — R K v|| M i ||v — xIImi, 

and thus ||v — RhvUh 1 ^ ||v— xIIh 1 for all x £ Sh- The first claim now follows 
by Theorem 14.11 

A duality argument [8j yields the second claim. Given v there is a 
unique u G H 1 satisfying (I + A)u = v — RhV, or equivalently (since A 
is self-adjoint) 

di (w, u) = (w, v — R h v) for all w G H 1 , 

Taking w = v — Rh.v and applying (14. 3p . we have for every x G Sh, 

(v - R h v, v - R h v) = ai (v - R h v, u) = m (v - R h v, u - x) 

< ||v - R h v|| h i||u-x||hi < Ch v_1 ||v|| H v||u-x||Hi- 

By Theorem 14.11 with q = 1 and y = 2 ^ 2t, there is a x G Sh such that 
1 1 t-L — xIIh 1 ^ Ch,||u|| H 2 , so 

||v — Rhv|| 2 ^ Ch, v ||v||H^||u|| H 2, 

and the result follows because ||u||h2 = ||(I + A)u|| = ||v — Rhv||. □ 
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4.2 Contour integral estimate 

We see from (jl,3J) and ( 13 . 1 9[) that, assuming Uoh = PhUo, 

u(z) = (zl + A) -1 g(z) and u H (z) = (zl + A H ) _1 Phg(z), 

so 

u H (z)-u(z) = G H (z)g(z) where G h {z) := (zl + A h ) _1 Ph — (zl + A) -1 . 

Deforming the integration contour in the Laplace inversion formula to V — 
9£p\ we can represent the error in the semidiscrete solution as follows: 



u H (t) -u(t) = ~~~T 
2m 



e zt G H (z)g(z)dz. (4.4) 



The next lemma allows us to estimate this integral. 
Lemma 4.4. // ^ v ^ 2t ; then 

||G H (z)v|| ^ Ch£||v|| H v-2, forzeL^ andveH"- 2 . 

Proof. Recall that £(z) := (zl + A)" 1 , and let £ h (z) := (zI + A H ) _1 . We 
split Gh(z) into two terms, 

Gh(z) = (P H -I)£(z)+ [£ h (z)Pk-Ph£(z)]. (4.5) 

Since A£(z) = (zl + A - zl)(zl + A)" 1 = I - z(zl + A) -1 , the resolvent 
estimate (12. ip shows that 

||A£(z)v|| < C||v|| forzGl^. 

Moreover, since (I + A)£(z) = I + (l — z)£(z) and since (1 + A) 1 ^ 2 commutes 
with (I + A)£(z), we have || (I + A)£(z)v|| H q ^ C|l — z||z| _1 ||v||h<i for any 
q £ K, and thus by Corollary 14.21 

||(P h -I)£(z)v|| < Ch v ||£(z)v|| H - = Cfi v ||(I + A)£(z)v|| H v-2 < Ch v ||v|| H -2, 

noting that |1 — z||z| _1 ^ Co,^ for zeI?. 

To estimate the second term in (14.51) . we write 

£ H (z)P h - Ph£(z) = £h(z)Ph(zI + A)£(z) - £ h (z)(zl + Ah)P h £(z) 
= £ h (z)[PhA-A h Ph]£(z). 
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For all u, w G H 1 , 



(P H (I + A)u,w) = ((I + A)u,P H w) = ai(u,P h w) = ai(R h u,P h w) 
= ((I + A h )RHU,P h w) = ((I + A H )RhU,w), 

so P H (I + A) = (I + A H )Rh and thus 

£h(z)Ph - Ph£(z) = £ h (z)(I + AJPh(Rh - I)£(z). 

Since £h(z)(I + Ah) = I + (1 — z)£h(z) the resolvent estimate (12.11) and 
Theorem 14.31 imply that 

||[£ h (z)P h -Ph£(z)]v|| < (1 + C|l - zllzr 1 ) || (R h - I)£(z)v|| 

< Ch v ||£(z)v|| H - = CH V ||(I + A)£(z)v|| H v- 2 < Ch v ||v|| H v-2, 

noting again that |1 — z||z| _1 < Cw,$ for z G lg\ □ 

Theorem 4.5. Let u be the solution of (II. ip anc? Zet Uh 6e £/ie semidiscrete 
approximation given by (I3.13P . // ^ v ^ 2x, then 

||u h (t)-u(t)|| ^Ch^e^dluollHv-a + HfllH-^), for t>0. 

Proof Let V± be the half-line z = cu + se ±l13 for < s < oo, so that 
T = r + — F_. Since £Rz = cu — cs where c = —cos (3 > 0, by applying 
Lemma [4.41 we have 



e G H (z)g(z) dz 



e (a, - cs)t ||G h (z)q(zl|| ds 



^ Ce^\\g\\ H , 



,— est 



ds, 



and the error bound follows at once from the integral representation (14. 4p . □ 

Combining Theorems 13.11 and 14. 5| we conclude that provided Uo and f 
have the appropriate spatial regularity, 

||U N)H (t) -u(t)|| = 0(lg(p r N)e^ N +H 2 X T ) for T/2 ^ t < 2T, (4.6) 

where the constant includes a factor (l+T _1 )e 2a)T . Moreover, in the next sec- 
tion (Theorem 14. 8| Part 2) we will see that when f = the error bound (14.61) 
remains valid even if the initial data is not regular. 
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4.3 Nonsmooth initial data 

Consider the case f = 0, that is, 

3 t u - A*u = on § n for t > 0, with u = u when t = 0, (4.7) 

and the corresponding semidiscrete problem in which : [0, oo) — > Sh 
satisfies 

3 t Uh — A^Uh = on § n for t > 0, with u = u h when t = 0, (4.8) 
where A^:Sh— Sh is defined by 

<- a h^)X) = a(ip,x) = (gradip,gradx) for all ^,xe S h ; 

compare with ( 13. 18ft . In contrast to the forgoing analysis, we now permit the 
initial data Uo to be an arbitrary function in L.2(§ TL ). 

By separating variables, we obtain an expansion in spherical harmonics, 

oo N(u,£] 

u(t) = £(t)uo = X L e- Att K) f , k Y £k , (4.9) 

1=0 k=l 

that implies the smoothing property in the next theorem. 

Theorem 4.6. Let < q < v and m G {0,1,2,...}. If u G H s then 
£(t)u G H v and 

||3™£(t)v|| H v < C T t- (v - q]/2 - m ||v|| H q, for < t ^ T. 

Proof. Adapting the argument of Thomee [T3| Lemma 3.2], we see from (14. 9 j) 
that the generalized Fourier coefficients of 9™£(t)uo are 

(ar£(t)uo,Y tv ) = (-A £ ) m e- AEt CM fk , 

so by (J23D, 

oo N(«,n) 

||3r£(t)u || 2 H . =]T_(l+\ t V\l m e- 2 ^ Y. I^Wf- 

t=0 k=l 

The result follows because, with s = A^t, 

f- q+2m (l + A £ r- q A 2m e- 2Aft < (T + s) v - q s 2m e- 2s < C T for < t ^ T. 

□ 
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Let 7 : L 2 — > H 2 be the solution operator for the elliptic problem 

u-A*u = f on § n , 

that is, Tf := u. Thus, 

Ql (Tf,v) = (f,v) forallveH 1 , 
and we can define Th : L 2 — > Sh by 

ai(T h f,x) = (f,x) for a11 X e S h . 
It follows that T h f = R H u = R h Tf and R H = 7 h {l - A*). Since 

(f , T h w) = ai (T hf , T h w) for all f , v e L 2 , 

we see that Th is self-adjoint and (taking w = f ) strictly positive-definite. 

Rewriting the homogeneous equation (14 .7|) as 3 t u + (I — A*)u = u, we 
see that 

T3 t u + u = Tu for t > 0, with u(0) = u , 
and similarly the corresponding semidiscrete problem (14. 8 P is equivalent to 

T H 3 t u h + u H = T h u h for t > 0, with u H (0) = Uoh- 

Thus, the error e = % — u satisfies 

T h 3 t e + e = T H e + p where p = (R H — I)u. (4-10) 

Lemma 4.7. With the notation above, if u h = Ph"U-o then 

||e(t)|| 2 < C T (j|p(t)|| 2 + -j- (s 2 ||3 s p|| 2 + ||p(s)|| 2 ) ds^j forO<t^T. 

Proof. We modify the argument of Thomee [13, Lemma 3.3]. Taking the 
inner product of (14.1 Oft with 3 t e gives 

(T H 3 t e, 3 t e) + (e, 3 t e) = (T h e + p, 3 t e), 

and since (T h 3 t e, 3 t e) ^ and (e, 3 t e) = (l/2)3 t ||e|| 2 , it follows that 

3 t ||e|| 2 ^ 2(T h e + p,3 t e), 
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implying that 

3 t (t||e|| 2 ) = ||e|| 2 + t3 t ||e|| 2 < ||ef + 2t(T h e + p, 3 t e). 

Since 



and 

we have 



2t(T h e,3 t e) = t3 t (T h e,e> ^ 3 t (t(T h e,e» 
t(p, 3 t e) = 3 t (t(p, e)) - t(3 t p, e) - (p, e), 



3 t (t||e|| 2 ) ^ ||e|| 2 + 3 t (t(T H e + 2p,e)) - 2t(3 t p, e) - 2(p, e), 
so integration gives 



tllelr ^ 



ft 



|e(s)|| 2 ds + t(7 h e + 2p,e) + 2 



I(s3 s p + p(s),e(s))| ds, 



and using 2(p,e) ^ 4||p|| 2 + (l/2)||e| 



t||e|| 2 ^ 2t(T H e,e)+8t||p|| 2 + 2 



t (s 2 ||3 s p|| 2 +||p(s)|| 2 + 2||e(s)|| 2 )ds. (4.11) 



To deal with the terms in e on the right-hand side, take the inner product 
of (I4.10p with e, obtaining 

U/2)3 t (T H e, e) + ||e|| 2 = (7 h e + p, e), 

or equivalently, 3 t (Th.e, e) — 2(Th.e, e) +2||e|| 2 = 2(p, e). After multiplying by 
the integrating factor e~ 2t , 



3 t (e- 2t (T h e,e)) +2e- 2t ||e|| 2 = 2e- 2t (p,e), 
and the choice u h = Ph u o means that Tne(O) = because 

(T h e(0), w) = (T H (P h - I)u , w) = ((P H - I)u , 7 h w) = 
for every w G L 2 . Thus, 



(4.12) 



rt 



e- 2t (7 h e,e) + 2 



ft 



e- 2s ||efs)|| 2 ds = 2 



e- 2s (p(s),e(s))ds 



Y 2s (||p(s)|| 2 + ||e(s)|| 2 ) ds, 
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implying that 



ft 



(T H e, e) 



,2(t-s)| 



efslll 2 ds < 



e 2 ^||p(s)|| 2 ds. 



Hence, 



ft 



2t(T H e, e) + 4 



|e(s)|| 2 ds 2max(t,2) 



t c 2 f^-3||p(fi)|| 2 ds, 



and inserting this bound in (14. lip gives 

2 



|e(t)|| 2 ^8||p(t)|| 2 + - 



s 2 ||3 s p|| 2 ds +3max(l,2t- 1 ] 



e 2(t - s) ||p(s)|| 2 ds. 



□ 



Theorem 4.8. Let u be the solution of the homogeneous problem ( 14. 7ft 
with initial data Uq, let % be the semidiscrete approximation given by (14. 8 p 
with u h = Ph u o- 1 ^ y ^ 2t: 

i. if uq e H v (S n ) ; &en 

||u H (t)-u(t)|| Ct^||u ||h^ /orO^t^T; 

<?. i/uo G L 2 (S> n ) and 2t is an integer, then 

||u H (t)-u(t)|| < C T t-/ 2 h-||u || forO<t<:T. 



Proof. We see at once from Lemma 14.71 that 

||e(t)||<C t sup (||p(s)||+s||a s p||), 

and if Uo G H v then, by Theorems 14.31 and 14.61 

||p(s)|| + s||3 8 p(s)|| <CH v (||u(s)|| H - + s||a s u(s)||m) < CIO 

which proves Part 1. 

Assume now that Uo € L 2 . By Theorems 14.31 and 14.61 



hp 



||p(t)|| = ||u(t)-R H u(t)|| < Ch||u(t)|| H i < Cht- 1 / 2 
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and the expansion (14. 9 j) in spherical harmonics implies that 



||p(s)fds < Ch 2 



|ufs)|| H i ds 



N(n,f ) 



Ch 2 Jjl+A £ ) Y_ |M 



Ck 



rt 



o 



-2A e s 



ds. 



£=0 k=l 

If £ ^ 1 then A^ ^ Ai = n so the substitution s = cr/A^ gives 



(1+A f 
and thus 



1 + A, 



rApt 



e - 2A ' s ds = 

o A 



e" 2a d0 < U + u- 1 ; 



e~ 2a dcr < 1, 



iip(s)ii 2 ds^ch 2 (ti(^r) 01 i 2 +^ y; iSwi 2 ) <c t h 2 Hii 2 . 

V £=1 k=l / 



Similarly, 



ft 



s 2 ||3 s p|| 2 ds ^ CH 2 



s ||3 s u(s)|| H i ds 



CH 2 ^(1+A £ )A 2 s 2 e- 2AtS ds \M 



fkl 



k=l 



and for all £ > 1, 



(1 + A,)A 2 



1 + A, 



s 2 e- 2AtS ds 
o A 



a 2 e 2a da^ C, 



so Jg s 2 ||3 s p|| 2 ds ^ Ch, 2 ||xio|| 2 . Applying Lemma I4T71 Part 2 follows in the 
special case v = 1. 

To deal with case v = 2t, we introduce the solution operator for the 
semidiscrete problem, £n(t)uo := Uh(t), and use the semigroup property: 
£(s + 1) = £(s)£(t) and £ H (s + 1) = £h(s)£h.(t) for all s and t. The error 
operator ^(t) = £n(t) — £(t) satisfies the identity 

5F H (t) - Jh(V2) 2 = £ H (t/2) 2 - £(t/2) 2 - [£ H (t/2) - £(t/2)] 2 
= 5F h (t/2)£(t/2) + £(t/2)5 h (t/2), 
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and by Part 1 and Theorem 14. 6[ 

||3 r h (t/2)£(t/2)u || < Ch v || £(t/2)uo || H - < Ch v (t/2)- v/2 ||u ||. 

Since £(t/2) and 3 r n(t/2) are self-adjoint in L 2 , the same estimate holds for 
the reversed product £(t/2)3 r H(t/2), and therefore 

H^hWuoII < Ct- v/2 h v ||u || + Ct- 1/2 K||3 r H (t/2)ix ||. (4.13) 

The stability estimates ||£(t)uo|| ^ ||xxo 1] an d ||£h.(t)uo|| ^ C||oxo|| mean that 
it suffices to consider the case t _1 ^ 2 h, ^ 1, when repeated application of the 
estimate (I4.13P gives 

||5 h (t)uo|| < Ct- v/2 h v ||u || + C(t- 1/2 h) j ||^H(t/2 j )u || 

for ) = 0, 1, 2, . . . , v = 2t, and thus H^hWixoII ^ Ct _T h 2T ||uo||. For the 
remaining case 1 < v < 2t, let = v/(2t) and observe that 

||5 h (t)uo|| = ll^hWuoll^HJHWuoll 9 

< C||u |r- e [(t- 1/2 ri) 2T ||uo||] e = Ct- Y/2 h Y ||u ||. 

□ 



5 Numerical experiments 

We present the results of some numerical experiments with two model prob- 
lems. In both cases, the integration contour ( 13. 4p and quadrature step k are 
chosen as in Theorem 13.11 with 

T = 1. cu = L 9 = 1/2, 5 = 7t/4, r = tt/4; 

Figure [T] shows the case N = 20. Our conference paper [3] presents some 
earlier numerical examples. 



5.1 A scalar problem 

Consider the ODE u' + u = f (t) for t > 0, with u(0) = 1. We choose the 
source term f so that the exact solution is 

4 t 3/2 

u(t) = 1 + — — 
3V7t 



20 




-80 1 1 1 1 1 1 1 

-60 -50 -40 -30 -20 -10 10 

Figure 1: The integration contour V and quadrature points Zj when N = 20. 

which has the Laplace transform u(z) = tT x +z~ 5 ^ 2 . In this case, no spatial 
discretization is required, and the numerical solution Un is given by ( 13. 8h . 
Table [1] shows the error at t = 2 for different values of N . The rapid conver- 
gence is consistent with the error bound of Theorem 13.11 but as N increases 
the quadrature eventually becomes unstable. 



N 


10 


20 


30 


35 


40 


|U N (2)-u(2)| 


1.71E-04 


6.44E-08 


3.75E-11 


7.52E-13 


1.16E-12 



Table 1: Errors for a scalar problem. 



5.2 Heat equation on the unit sphere 

Fix < a < 1 and define u : § 2 — > C for x = (x 1; x 2 , x 3 ) G § 2 by 

Jl, if a< x 3 < 1, 
UoM = <^ (5.1) 
10, if — 1 < x 3 < a. 
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m 


PmM 


Smoothness 


T 


2 


(l-r)^(3 + 18r + 35r 2 ) 


C 4 


7/2 


3 


(1 - r)%{l + 8r + 25r 2 + 32r 3 ) 


c 6 


9/2 



Table 2: The compactly supported SRBFs of Wendland [16J. 



This axially symmetric function has the Fourier-Legendre expansion 



Uo(x) = Y_ ( u o)cJM*3), where [u ) t = 



21 + 1 



£=0 



Ppftldt. 



The zeroth coefficient is (uo) = (1 — a)/2, and the remaining coefficients are 
expressible in terms of Jacobi polynomials [1, page 172], [7, Formula 18.9.15], 

Q < = HT^ eIFTT) p;(a) = (1 - ^^W 1 p! " 1(a) for ( > 1; 

consequently (u^) e = 0(£~ 1/2 ) as I -»■ oo [12, Theorem 7.32.2]. 

The PDE u t — A*u = with initial data ( 15. ID describes heat diffu- 
sion from a spherical cap about the north pole onto the surface of the unit 
sphere S 2 . By separating variables, we find that the exact solution is 

oo 

u(x,t) = Y_ e~' (<+1)t (uo) e P«(x 3 ), for x = (x 1 ,x 2 ,x 3 ) G § 2 . 

£=0 

For the spatial discretization, we use the compactly supported radial basis 
functions introduced by Wendland [16], for which the strictly positive-definite 
kernel has the form 



®{x,y) = p m (v/2-2x-y). 

In Table |2j we show p 2 and p3 explicitly, along with the values of the expo- 
nent t in (12. 81) . We generate the set of points X using an equal area parti- 
tioning algorithm of Saff and Kuijlaars [9]. To compute the inner products 
arising in the matrix entries (13.151) and the load vector components G p (z), 
we use a quadrature approximation of the form 

2tt R R/2 

v dS « — 2~_ /_ w p v (si n 9p cos 4> q , sin 9 p sin (b q , cos 9 p ) , (5.2) 
R — ■ — ■ 

q=lp=l 
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K 


200 


400 


600 


801 


1001 






H x 


0.1796 


0.1281 


0.1039 


0.0888 


0.0794 






R 


200 


200 


200 


500 


500 


N 


= 10 


6max 


5.67E-05 


5.06E-06 


1.67E-06 


9.89E-07 


8.24E-07 






e 2 


4.63E-05 


3.60E-06 


1.10E-06 


8.23E-07 


7.81E-07 






EOC(e 2 ) 




7.56E+00 


5.66E+00 


1.84E+00 


4.66E-01 


A. 1 

N 


= 20 


6max 


5.61F-05 


4.47h-06 


1.03b-06 


3.26F-07 


1 AO T~? r\ T 

1.48F-07 






e 2 


4.61E-05 


3.53E-06 


8.07E-07 


2.64E-07 


1.20E-07 






EOC(e 2 ) 




7.60E+00 


7.05E+00 


7.11E+00 


7.06E+00 


N 


= 30 


6max 


5.61E-05 


4.47E-06 


1.03E-06 


3.26E-07 


1.48E-07 






e 2 


4.61E-05 


3.53E-06 


8.07E-07 


2.64E-07 


1.20E-07 






EOC(e 2 ) 




7.60E+00 


7.05E+00 


7.11E+00 


7.06E+00 


N 


= 35 


6max 


5.61E-05 


4.47E-06 


1.03E-06 


3.26E-07 


1.48E-07 






e 2 


4.61E-05 


3.53E-06 


8.07E-07 


2.64E-07 


1.20E-07 






EOC(e 2 ) 




7.60E+00 


7.05E+00 


7.11E+00 


7.06E+00 



Table 3: Numerical results with SRBFs constructed using p 2 . 



for an even number R ^ 2, where J_ 1 f (z) dz ~ 2^p= 2 i w p f (cos 9 p ) is a Gauss- 
Legendre rule and cb q = 27iq/R. The error in the approximation f)5.2p is zero 
if the integrand v is a polynomial of total degree R — 1 or less. 

In the numerical experiments, we let a = 0.9 in the definition ( 15.11) of Uo- 
Tables [3] and H] show values of the quantities 

e ma x = max|U NiH (x, 1) - u(x, 1)1 



( Y_ w x |U Njh (x, 1) - u(x, 1) 



1/2 



and 

e 2 

SceQ 

for different choices of K and R. Here, Q is the set of quadrature points. 

Since Uo € L 2 (§ 2 ), we expect from Theorem 14.81 and the triangle inequal- 
ity fl3~2H that if N is sufficiently large then e 2 = 0(h 2T ) — that is, 0(h 7 ) 
using p 2 , and 0(fi 9 ) using p 3 . The observed convergence rates are close to 
these predicted values. We remark that when K = 1001, the condition num- 
ber of the linear system (13 . 1 9[) is around 10 7 using p 2 , and around 10 9 using 
P3, so we cannot expect to reduce the error much below the smallest values 
shown in the tables. 
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K 


200 


400 


600 


801 


1001 




H x 


0.1796 


0.1281 


0.1039 


0.0888 


0.0794 




R 


200 


200 


200 


500 


500 


N = 10 


6max 


6.86E-05 


3.84E-06 


1.17E-06 


8.35E-07 


7.79E-07 




e 2 


3.93E-05 


1.63E-06 


7.85E-07 


7.74E-07 


7.71E-07 




EOC(e 2 ) 




9.41E+00 


3.50E+00 


9.13E-02 


3.37E-02 


N = 20 


6max 


6.78L-05 


3.11L-06 


4.54L-07 


8.98L-08 


O O A T? n o 

3.24L-08 




e 2 


3.91E-05 


1.45E-06 


2.12E-07 


4.83E-08 


1.73E-08 




EOC(e 2 ) 




9.75E+00 


9.17E+00 


9.41E+00 


9.18E+00 


N =30 


6max 


6.78E-05 


3.11E-06 


4.54E-07 


8.98E-08 


3.24E-08 




e 2 


3.91E-05 


1.45E-06 


2.12E-07 


4.83E-08 


1.73E-08 




EOC(e 2 ) 




9.75E+00 


9.17E+00 


9.41E+00 


9.18E+00 


N =35 


6max 


6.78E-05 


3.11E-06 


4.54E-07 


8.98E-08 


3.24E-08 




e 2 


3.91E-05 


1.45E-06 


2.12E-07 


4.83E-08 


1.73E-08 




EOC(e 2 ) 




9.75E+00 


9.17E+00 


9.41E+00 


9.18E+00 



Table 4: Numerical results with SRBFs constructed using p 3 . 
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